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Abstract. We study the effects of a probabilistic refractory period in the collective 
behavior of coupled discrete-time excitable cells (SIRS-like cellular automata). Using 
mean-field analysis and simulations, we show that a synchronized phase with stable 
collective oscillations exists even with non-deterministic refractory periods. Moreover, 
further increasing the coupling strength leads to a reentrant transition, where the 
synchronized phase loses stability. In an intermediate regime, we also observe 
bistability (and consequently hysteresis) between a synchronized phase and an active 
but incoherent phase without oscillations. The onset of the oscillations appears in 
the mean-field equations as a Neimark-Sacker bifurcation, the nature of which (i.e. 
super- or subcritical) is determined by the first Lyapunov coefficient. This allows us 
to determine the borders of the oscillating and of the bistable regions. The mean-field 
prediction thus obtained agrees quantitatively with simulations of complete graphs and, 
for random graphs, qualitatively predicts the overall structure of the phase diagram. 
The latter can be obtained from simulations by defining an order parameter q suited 
for detecting collective oscillations of excitable elements. We briefly review other 
commonly used order parameters and show (via data collapse) that q satisfies the 
expected finite size scaling relations. 
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1. Introduction 



Understanding collective oscillations of coupled nonlinear elements remains a challenge 
from both theoretical and experimental viewpoints. From the theoretical side, much 
progress has been accomplished since the seminal works of Winfree and Kuramoto [H 
El [31 m [5] on coupled oscillators, upon which recent literature has expanded to include 
effects of e.g. complex topologies |E] and noise [71 El [9l HOI |TT1 |T2l [13]. From the 
experimental side, the subject has a longstanding importance in neuroscience: collective 
neuronal oscillations stood for a long time as candidates for a solution of the so-called 
binding problem [14], but emphasis has recently shifted to attentional processes [T5] . 

Here we are interested in collective oscillations of units which are excitable, i.e. 
not intrinsically oscillatory. This topic has been experimentally observed in a variety 
of scenarios, from neuroscience [161 E] to chemistry [H], but theoretical approaches 
have been relatively scarce [TH [201 Ell [221 [231 EH [25]. In particular, it is not entirely 
clear to which extent these oscillations are robust with respect to noise. On the one 
hand, recent studies have shown sufficient conditions for the onset of global oscillations 
of deterministic excitable units with noisy coupling, emphasizing e.g. the interplay 
between coupling strength and characteristic time scales of the units [201 EI], or the 
importance of the topology of the network [19]. On the other hand, the susceptible- 
infected-recovered-susceptible (SIRS) model on a lattice (i.e. a minimum three-state 
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excitable model) has been thoroughly studied, with several results strongly suggesting 
that its stochastic Markovian version does not yield sustained oscillations [261 EH [281 123] 
(for non-Markovian models, see e.g. [291 130] ). 

This raises the question whether it is possible to find sustained global oscillations 
in a network of excitable units whose intrinsic dynamics (not only the coupling) is 
non-deterministic. We therefore propose and study a simple probabilistic model which 
has a well-defined deterministic limit. To study the phase transitions in the model, 
we employ an order parameter specifically tailored to assess collective oscillations of 
excitable systems. These are described in section [2l We study complete graphs as 
well as random graphs, comparing mean-field calculations with simulations. Results are 
presented in sections [3] and HI while section [5] brings our concluding remarks. 

2. Model 

2.1. Excitable cellular automata 

The minimum model of an excitable system consists of three states, representing 
quiescence (state 0), excitation (state 1) and refractoriness (state 2) [29] (a prototypical 
example being the SIRS model). As shown by Girvan et al., however, such a cyclic three- 
state deterministic cellular automaton fails to exhibit sustained collective oscillations. 
For them to become stable, their model needs at least 2 refractory states [20j. To obtain 
an arbitrary number of refractory states, for each site j = 1, N, let Sj = 0, 1, 2, r be 
the consecutive states of the unit (out of the r + 1 states, the last r — 1 are refractory [20] . 
see figure H]). 

The cellular automaton version of the probabilistic SIRS model (also called the 
probabilistic Greenberg-Hastings model pT]) corresponds to r = 2, with intrinsic 
transitions 1 — )■ 2 and 2 — > governed by constant probabilities, whereas the 
transition — )■ 1 occurs with a probability that usually increases linearly with the 
number of excited neighbors [32j (for a study with nonlinear coupling, see [23]). 
In the model studied by Girvan et al., on the other hand, all intrinsic transitions 
{sj — 7- {sj + 1) mod (r + 1), sj ^ 0) are deterministic. 

Here we study an intermediate variant of these models, where all intrinsic transitions 
are deterministic but the last one, which occurs with probability (see figure H]). The 
idea is to have a minimum model (lest the number of additional parameters becomes 
too large) which incorporates non-determinism in the intrinsic dynamics. The choice to 
make the transition from the last refractory state probabilistic is natural and comes from 
neuroscience: neuronal dynamics depend on ionic channels which are stochastic [331 IM]. 
so that a neuron may or may not fire when stimulated at the end of its refractory period 
(the so-called relative refractory period) [33] . 
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Figure 1. Single-cell dynamics, pinf is the probability of activation (0 1) from 
neighbours, is the probability of transitioning from the relative refractory state (r) 
to the rest state (0). Light gray states are refractory. All other transitions are 
deterministic. 



2.2. Coupling 

The only transition that still needs to be described is the excitation process — )■ 1. We 
assume each site j is symmetrically connected with kj other sites. Each active site has a 
probability a/K of activating a resting neighbour, where cr is a control parameter (which 
corresponds to the system branching ratio [35l [22]) and K is the average connectivity 
{K = {kj)). The fraction of active sites 

1 ^ 

is used to measure network activity at time t. At cr = 1 the model shows a transition 
from an absorbing to an active state, but without sustained oscillations [22j. With 
deterministic units {p-y = 1), the system undergoes a transition to the oscillatory regime 
at 0" = adP'y = 1) > 1, which persists indefinitely if a is further increased pOl 122] . 

To motivate the analysis to be developed in section [HI figure [2] shows examples 
of single-run results in a complete graph with = 5 x 10^ for p^ = 0.85 and 
increasing values of cr. For probabilistic units, we observe the transition to an active but 
nonoscillating state at cr = 1 [Figure [2]^a)-(b)] and the second transition to an oscillating 
state for larger a [Figure [2]^b)-(c)]. Contrary to what is observed in the deterministic 
model, however, in the probabilistic model this second transition is reentrant with 
respect to the coupling strength cr, as shown in figure [2]^c)-(d). 

The nature of this reentrant transition will be clarified in section [3] In order to 
analyze it properly, though, one needs to define an adequate order parameter to detect 
synchronization among excitable elements. 



2.3. Order parameters 

Most studies of synchronization employ the Kuramoto order parameter [21 [3], which 
corresponds to the time and ensemble average of the norm of the complex vector 

1 ^ 

Z{t)^j^J2^'''^'^ ^ (2) 

where 6j = 2TTSj/ {t + 1). Note that Z corresponds to the center of mass of the phases of 
the units. This works fine when the system is composed of coupled uniform oscillators. 
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Figure 2. Time series for — 0.85 and increasing values of a on a complete graph 
with = 5 X 10'^ sites. Respectively for increasing a (top-bottom): (a) absorbing 
(not-active) state, (b) active state without oscillations, (c) with oscillations and (d) 
again active without oscillations. Po(0) — 0.95, /o(l) = 0.05. r = 3. 



because rotational symmetry will ensure that, in the absence of sustained collective 
oscillations, the order parameter vanishes in the thermodynamic limit. Consider, 
however, the present case of excitable elements. The trivial absorbing state {sj = 0, 
Vj), which is always a collective solution of the dynamics, yields a nonzero (in fact, 
maximum!) Kuramoto order parameter. Indeed, in the absorbing state units are 
all "perfectly synchronized" in the sense that they always have the same state. But 
this is clearly not what one wants to detect. This problem persists for cr > 1 (below 
the onset of collective oscillations), where a small fraction of the units are active (on 
average), whereas most remain quiescent. In that case, Z has a constant bias towards 
the absorbing state, around which it will fluctuate. 

Collective oscillations correspond to rotations of Z which may be misdetected in 
the averaging procedure owing to a lurking constant vector. One possibility which has 
been used to avoid the weight of the absorbing state is excluding terms with Sj = 
from the sum in eq. (j2]) [191 El] • Another strategy makes use of the standard deviation 
(measured along time) of -Pt(l) to detect oscillations [36j, though neither procedure is 
easily extensible to continuous-phase systems. 

To account for a system of continuous-phase units which may have an arbitrary 
number of preferred phases, one could employ the angular momentum L = XdtY — 
YdtX, where Z = X + iY [371 EH]- In our cellular automata, this would require the 
discretization of the time derivative, which could introduce unnecessary numerical errors. 
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(3) 



which amounts to subtracting the constant bias from Z before the averaging 
procedure [39]. However, this can be computationally expensive, requiring the storage 
of the whole time series. Making use of a similar idea, but with much less computational 
bookkeeping, in the following we will characterize collective oscillations via the order 
parameter 



which can be seen as a generalized standard deviation of Z{t). Differently from g, 
obtaining q is computationally inexpensive, as the means over time are now separated 
and may be calculated along with the simulation. In section 14.31 we will show that q 
satisfies scaling relations near a phase transition, as expected for an order parameter. 

One may grasp intuition about q by considering the different time series in figure [2j 
Let Pi be the stationary value limt^oo(-Pt(l))t, which is an order parameter in its own 
right, measuring whether or not the network is active [22]. In figs.[2](a) and (b), although 
the system goes from P* = to P^* 7^ as the system goes from an absorbing to an 
active fixed point, both have g = 0, because there are no oscillations after the transient. 
In figure [2]J^c), on the other hand, we have both P^ 7^ and g > 0, as the oscillatory state 
becomes stable after increasing a. Finally, for figure [2]^d), the oscillations are unstable, 
PI ^ and g = 0. 

3. Complete graph 

We start with the complete graph because it is presumably the topology most prone 
to exhibiting stable collective oscillations. Besides, it allows a comparison between 
simulations and analytical results (see below). From now on, we will focus on the effect 
of the probabilistic dynamics (p-y) on synchronization and will fix r = 3. 

3.1. Simulations 

We have simulated complete graphs {kj = K = N — 1) with sizes varying from = 10^ 
to = 10^. At each time step, intrinsic transitions occur as described in section O The 
transition — )■ 1 is governed by the number A'^i(t) of active (s, = 1) sites at time t, each 
of which can activate a quiescent cell with probability a/N. Therefore, the probability 
of a quiescent cell being activated by at least one of its Ni{t) active neighbours is 



which renders the simulations relatively simple despite the large system sizes. Finally, 
initial conditions must be chosen to avoid large amplitudes during the transient, which 
could throw the system into the absorbing state [20]. Apart from that, the effects we 





(4) 




(5) 
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report below are robust with respect to the initial conditions and we have arbitrarily 
fixed Po(0) = 0.8 and Po(l) = 0.2. 

To illustrate the kind of reentrant transition exemplified in figure El we show in 
figure [3]^a) the order parameter g as a function of the coupling parameter a. Starting at 
some amin > 1, for each value of a we let the system evolve during a transient of Urans 
time steps, after which we start measuring the order parameter g up to t = t^^ax time 
steps. We then increase cr by a constant amount 5a and repeat the procedure, with the 
initial condition of the system for each value corresponding to the final condition of the 
preceding value. Both constants are chosen so the rate of change is small {5a/tmax ^ !)• 
After a maximum value amax is reached, a is sequentially decreased by the same amount 
5a down to amin- 

As shown in figure [3|a), the reentrance of the transition to collective oscillations 
is captured by the non-monotonic behavior of q{a), which departs from zero at some 
lower value acij)-^) and returns to zero at some upper value a'^iP'-^)- Moreover, we have 
found that while the first transition is always continuous, the second transition can be 
discontinuous. The fingerprint of the discontinuity is the hysteresis observed in the order 
parameter: above a 2-, the only stable state of the system has constant (but nonzero) 
Pt(l), thus no oscillations (g = 0). If we decrease a, oscillations do not reappear at cig, 
but rather at a lower value a\. There is therefore a region of bistability a G ['7i,o"2] 
in parameter space where collective oscillations (g > 0) can coexist with an active 
(Pt(l) > 0) but non-oscillating (g = 0) state. As it turns out [see Figure |3]|^a)], for 
= 0.95 the size of this bistable region is rather sensitive to the system size. Smaller 
systems tend to be perturbed away from collective oscillations by larger fiuctuations, 
leading to smaller hysteresis cycles. 

As is decreased, the mean and variance of the refractory periods of the units 
increase, rendering the whole system noisier. This might be the explanation for the 
result in figure [3]^b), which shows a smaller reentrant region with oscillations (g > 0). 
The width of the hysteresis cycle also decreases, with both a^ and decreasing. 
Furthermore, (along with the width of the hysteresis cycle) becomes less sensitive 
to the system size, which could be due to the variance of the refractory periods 
overcoming the effects of small-size fluctuations. Albeit subtly, cxc also increases slowly 
with decreasing will be seen in flgure HI 

Further decreasing p^ [Figure [3]^c)], the bistable region vanishes, whereas the 
transition to a collectively oscillating state remains. Finally, for sufficiently small p^, 
collective oscillations are no longer stable [Figure |3]^d)]. 

In the following, we will show that these transitions can be quantitatively 
reproduced by a low- dimensional mean-fleld analysis. For a controlled comparison, 
(Tc and a^ were heuristically deflned as the values of a (averaged over n runs) where 
g first rose above some threshold value gmm oc I/a/ZV, respectively for increasing and 
decreasing values of a. On the other hand, was defined as the value of (increasing) 
a where g first fell below qmin- 
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Figure 3. Hysteresis loop for the mean-field solution (to be presented in section [ 
and complete graph {N = 10^ 2.5 x 10^ and 10^) with = 0.95, 0.9, 0.8 and 0.75. 
tmax = 1 X 10'^ steps {ttrans = 500) and 5a = 0.05. Mean (symbols) and standard errors 
(bars) calculated over 10 runs. Insets zoom into the interesting region 4 < cr < 10. 
Note that fluctuations in the inset of (d) decrease with increasing system size, where 
the mean-field approximation f section 13. 2p predicts 5 = in the thermodynamic limit 
(see also section ITS]) . 



3.2. Mean- field analysis 

In the following we apply the standard mean-field (MF) approximation [ID] to the 
equations governing our system. We follow closely the steps of Refs. [20, SUES], where 
every site is considered to have K neighbors, a fraction of which is excited at time 
t. If a given site is at rest [with probability -Pt(O)], the probability of it becoming excited 
by at least one of its excited neighbors is 

p.„,M=l-(l-!^)^ (6) 
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The dynamics of the system is then described by the following closed set of equations: 
Pi+i(0) = p,P,{t) + (1 - P,„;(t)) P,(0) (7) 
Pt+i(l) = P,„/(t)Pj(0) (8) 



(9) 

(10) 



Pt+^{s) =Pt{s-l) (2<s<r-l) 
Pt+i(r) = Pi(r-l) + (l-p^)Pt(r), 
where the normahzation condition 

r 
s=l 

renders ([7]) redundant and reduces the system to a r-dimensional map [H]. Therefore, 
increasing the duration of the refractory period amounts to an increase in the complexity 
of the mean-field calculations. 

For the complete graph, K = N — 1 and mean field is exact. In the thermodynamic 
limit, (El) becomes 



inf 



lim P 

Af-*-oo 

which, from (iTl) and ffTTl). leads to 



P 



t+i 



1) = (i-e-'^^'^i): 



s=l 



(12) 



(13) 



In other words, in the mean field approach the state of the system is completely described 
by Pt = (Pt(l), Pt(2), . . . , Pt(r))^, which evolves according to Pt+i = F{Pt). Note that 
Pi is the only component of F which is nonlinear [see (fT3l) ]. 

From ([9]), the fixed point for any 2<s<r — lis clearly P* = Poo(s) = Poo(s — 1) = 
■ ■ ■ = P* which, upon substitution in ffTOl) . gives P* = Pf/p^. Finally, in its steady state, 
(fT3l) becomes: 



P* 



:i 



(14) 



- 1 + — P* 
PiJ . 

which can be numerically solved. Expanding near Pj* ~ (note that Pj* = is always 
a solution), one easily obtains the transition from an absorbing to an active (steady) 
state at 0" = 1 [22j. The collective oscillations appear when the active state becomes 
unstable. 



3.3. Linear stability 

Considering a small perturbation rit{s) such that Pt(s) = P* + 77t(s), the linearized 
dynamics can be written as rff+i = Afft, where Aij = dFi/dPt{j)\p* is the Jacobian 
matrix calculated at the fixed point: 



/ gia,P^) (e-'^^i -1 
1 
A= 1 

V 



T _ 1 



-<7P* 



1) \ 













(15) 
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and gicr,Pi) = e"^^* - 1 + ae-'^-^i [1 - (r - 1 + l/p^)P*]. The eigenvalues {/ij}J=i 
of A determine whether P* is stable (maxj < 1) or unstable (maxj > 1) (for 
simplicity, in the following we employ /i = /x^, where k = argmaxj|/ij|). 

We expect to pinpoint the transition to the oscillatory state by looking for a 
Neimark-Sacker (NS) bifurcation in the mean field equations (which is the discrete- 
time analog of the Andronov-Hopf (AH) bifurcation in continuous time [12] )• In other 
words, for fixed p^, we have = 1 with Im(/i) 7^ at cr = a^^. The relation between 
a^^ and the pair {aj, cTg} of critical values depicted in figure [3]^a) will be clarified below. 

Like the AH bifurcation, the NS bifurcation also comes in two different flavours: in 
the supercritical stable closed invariant curve (CIC — the discrete-time analog 

of a limit cycle) is born at a^^ and grows continually from zero amplitude; in the 
subcritical case, an unstable CIC exists below cr^'^ and engulfs the fixed point P* at 
(7^^ (above which the system is typically attracted to another pre-existing, but stable, 
CIC). Since the order parameter q increases with the amplitude of the oscillations (which 
is, roughly speaking, proportional to ||r/||), a supercritical (subcritical) NS bifurcation in 
the mean-field equations is suggestive of a continuous (discontinuous) phase transition 
in the system (see e.g. [23]). 

The sign of the first Lyapunov coefficient h [12] indicates if the Neimark-Sacker 
bifurcation is supercritical (/i < 0) or subcritical [li > 0). Its calculation is briefly 
reviewed in the Appendix. We now have the necessary tools to unveil the complete 
phase diagram. 

3.4- Phase diagram 

We have run simulations of complete graphs with = 10^ excitable units and employed 
the protocol described in section [3A] with tmax = 10^ to detect the width of the hysteresis 
loop (coexistence region). We have tested and verified that longer values of tmax do not 
change our results significantly. The phase diagram thus obtained from the simulations 
is shown with symbols in figure H] (the horizontal gray lines show the values of used 
in figure E]). To obtain the NS lines of the mean- field equations, we have numerically 
explored a special test function |32] ^Nsi'^iP-r) = Y[m<n(^ ~ Z^n/^m), which changes its 
sign at the NS bifurcation. As the first Lyapunov coefficient /i determines whether the 
bifurcation is super- or subcritical, its value along the bifurcation line is shown in the 
inset of figure H] (changing sign at ctt)- The solid lines in figure H] show the supercritical 
(red) and subcritical (blue) bifurcation curves where ^Nsi^r'^^'iP-y) = 0. 

Comparing the solid lines with the symbols in figure HI we observe that linear 
stability analysis accurately accounts for the transition from an active (but non- 
oscillating) phase to an oscillating phase. In other words, it correctly predicts the 
lines (JciP'y) and al{p^), where in both cases the non-oscillating phase loses stability. 

The transition at o'2{p-y), however, cannot be predicted by linear analysis. Note 
that in this case it is the stable CIC that loses its stability (that of the fixed point 
P* remaining intact). This hints at the existence of a global bifurcation, which can be 
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Figure 4. Phase diagram for the complete graph (triangles) with N — 10^ (mean 
over 50 runs, t,nax — 10'^, Urans = 5 X 10^). Solid red (blue) line is the supercritical 
(subcritical) Neimark-S acker bifurcation predicted from linear analysis. Black line 
marks the (discontinuous) stability limit of the oscillating phase (as predicted by mean 
field). Note that the black and blue lines approach each other very closely before 
merging at ctt. Inset: first Lyapunov coefficient ?i((t). 



numerically detected in the MF equations by direct iteration of the map determined 
by equations ([9]), ( ITOl) and ( fT3|) . To compare the r-dimensional MF map with system 
simulations, we rewrite the complex vector Z from (|2]) and (fTTj) as 



where 0s = 2tts/{t + 1). We can thus calculate q for the MF map and subject it to the 
same protocol used for detecting the coexistence region in the simulations. The black 
solid line in figure H] shows the (J2{p-y) obtained by iteration of the map, which is in good 
agreement with simulations (symbols). 

Note that in the lower part of the oscillating phase {pj < 0.8) the order parameter 
detects oscillations in the simulations which are not predicted by the mean-field analysis. 
This phenomenon is due to stochastic oscillations, as recently explained by Risau- 
Gusman and Abramson [13]: the fixed point in the conflicting region is in fact stable, 
but with an eigenvalue with a nonzero imaginary part. Inevitable fluctuations throw 
the system away from the stable point, to which it returns in spiral-like trajectories, 
yielding a nonzero q even for very large system sizes [lH [23] . 

For = 1, we recover a quenched variant of the model by Girvan et al. [20j. In 
this regime where intrinsic transitions are deterministic, increasing the coupling will 



r 




(16) 
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only reinforce collective oscillations, and the fixed point P* never regains stability (i.e. 
al —7- oo). This suggests that even small amounts of noise in the intrinsic dynamics of 
excitable elements can lead to qualitatively different collective behavior in a regime of 
strong coupling. 

4. Random graph 

4-1. Mean- field and simulation results 

To understand network topology effects on synchronization, we study a bidirectional 
random graph similar to an Erdos-Renyi [H] network, where NK/2 links connect 
randomly chosen pairs p2j and remain frozen ("quenched") throughout each run (and 
in each run, a new realization of the network is created). The main difference with 
respect to the complete graph (CG) is that in the random graph (KG) the value of a 
is bounded from above: a < K [see ([6])]. The mean-field calculations, however, are 
otherwise similar to that of the complete graph, with ([6]) replacing ( fT2l) . We therefore 
applied to the RG problem the same procedures for determining the stability of the 
solutions, the nature of the NS bifurcation and the boundary of the bistability region 
(see section [3]). 

Given their uncorrelated assigment of links and short distances among sites, random 
graphs are usually regarded as the natural topology in which mean-field predictions are 
expected to hold. Indeed, simulations and mean-field calculations agree nearly perfectly 
as far as the phase transition at cr = 1 [22] is concerned. In figure [5]J^a)-(b) we see that 
a good agreement is also observed in the transitions to a synchronized phase for large 
values of K. Note, however, that simulations and mean-field predictions differ at the 
rightmost boundary of the bistable region, and the disagreement worsens as K decreases. 
As shown in figure [5]J^c), for smaller values of K bistability was not even detected in the 
simulations, and the oscillating phase is substantially smaller than predicted by mean 
field. 

4-2. Annealed random graphs 

Could correlations (which are neglected by the mean-field approximation) account for 
the discrepancy observed in figure [5]? In order to assess the role of the correlations 
associated with the quenched connectivity, we studied an annealed variant of the model 
where the K neighbors of each site are randomly chosen at each time step [20l HS]. 
Results are shown in figure [6l in which we restrict ourselves to smaller values of K 
because results are essentially indistinguishable from the quenched case for large K. 

For smaller three features are noteworthy. First, the agreement with mean 
field results is recovered (apart from the stochastic oscillations in the lower end of the 
oscillating phase, like in the previous cases). This therefore confirms the suspicion that 
correlations associated to the quenched connectivity can indeed undermine collective 
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Figure 5. Phase diagram for (quenclied) random graplrs witlr (a) K — 300, (b) 
K — 150 and (c) K — 30. Solid red (blue) line is the supercritical (subcritical) 
Neimark-Sacker bifurcation predicted from linear analysis. Black line marks the 
(discontinuous) stability limit of the oscillating phase (as predicted by mean field). 
Symbols are obtained from simulations (mean over 5 runs) with tmax = 3 x 10"^ steps 
{ttrans — 2 X 10^ stcps) and N — 10^. Standard errors are smaller than symbol size. 
The purple dashed line is a guide for the eyes and marks the stability limit of the 
oscillating phase for simulations. 



oscillations. This is not surprising, given the difficulty of establishing collective 
oscillations of excitable elements in hypercubic lattices [T9| [23]. 

Second, the bound a < K impoverishes the repertoire of detected phenomena for 
small K [see the forbidden gray regions in figure [6]J^b) and (c)]. Note that the coexistence 
region shrinks as K decreases. In fact, since varies very slowly with K [see the inset 
of figure El^c)], there is a minimum value of K [satisfying = ariKc)] below which 
the system shows no bistability (i.e. there is no change in the sign of li). We have 
numerically estimated Kc = 7.8074(1). 

Finally, note that the oscillating phase extends into lower values of for decreasing 
K. This is a rather counter-intuitive result. It means that, for fixed a and p^, it is 
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Figure 6. Phase diagram for (annealed) random graplrs witlr (a) K — 30, (b) K — 2Q 
and (c) K — 10. Solid red (blue) line is the supercritical (subcritical) Neimark- 
Sacker bifurcation predicted from linear analysis. Black line marks the (discontinuous) 
stability limit of the oscillating phase. Symbols are obtained from simulations (mean 
over 5 runs) with t^ax = 3 x 10'^ steps [ttrans — 2x 10^ steps) and N = 10^. Standard 
errors are smaller than symbol size. Grey shaded areas correspond to a forbidden 
region where a > K. Inset: ctt for different values of K, showing the existence of a 
Kc = ariKc) = 7.8074(1) below which no bistable region exists. 



possible to take the system from a non-oscillating to an oscillating phase by lowering 
the connectivity K. This is a particularity of the annealed RG (and the mean-field 
approximation), however. Note in figure [5] that the opposite (and expected) trend is 
observed for quenched random graphs. It remains to be studied whether refining the 
approximation (including e.g. first-neighbor correlations [HI [28]) can reconcile mean- 
field with quenched random graph results. 

4-3. Finite-size scaling 

Our phase diagrams rely heavily on the proposed order parameter q. The inset of figureO 
indicates that near cTc its behavior becomes increasingly abrupt with increasing system 
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size A^. In order to confirm that q indeed possesses the basic properties of a bona fide 
order parameter, here we show that at the supercritical NS bifurcation (i.e. a second 
order phase transition) it satisfies the scaling relations that would be expected from 
standard finite-size scaling (FSS) theory. 

Defining A = a - cTc, FSS predicts that q oc L^^/^^f {AL^^"^) for a lattice with 
linear size L [30], where the critical exponents are defined as q oc jAj'^ and C, oc |A|~'^-l 
in the limit N = L'^ ^ oo (where ^ is the correlation length). This holds for d below 
the upper critical dimension dc- For d > dc, mean-field exponents are expected and the 
scaling relation has to be modified [16] with L — )■ N'^^'^", so 

This modified version of the usual FSS relation was shown to also hold for infinitely 
coordinated networks |3Z] (i.e. the complete graph). 




N=7!00x10^ 
N^1, 0,0x1 0'^ 



10° 10^ 10^ 10^ 

Figure 7. Collapse in a probabilistic (p^ = 0.9) regime for K = 150, tXc = 5.16. We 
used /3 = 1/2, i'± = 1/2, dc = 4. Mean over 15 runs, with standard errors smaller than 
symbol size. Other parameters were tmax = ^0^,ttrans = 7 x lO'^, Sa ^ 5 x 10~^ 



For large argument, the scaling function in Equation f|T7|) becomes f{x) oc , 
as usual [30]. In the subcritical regime (A < 0), one expects q ~ 0{N-^/^) [3], 
so /{AN^/'^'^"^) oc A^-i/2jv/3/rfc;^x_ this to be true, f{x) oc x^"^ when x < 0. 
Figure [7] shows an excellent data collapse for (quenched) random graphs with different 
system sizes. Consistent power laws are obtained with standard mean-field exponents 
(as expected for random graphs), namely /3 = l/2[3],i/^ = l/2 and dc = 4: [30]. Similar 
results are obtained with complete graphs (not shown). 
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5. Concluding remarks 

We have studied the effects of a probabilistic refractory period in the collective behavior 
of a large number of coupled excitable cellular automata. We have obtained the mean- 
field solution of the model and compared it with simulations of complete as well as 
random graphs. The continuous phase transition to a synchronized regime is associated 
to a Neimark-S acker bifurcation in the mean-field equations. 

This scenario is similar to what has been previously obtained by Girvan et al. [20] 
in a model of deterministic excitable automata (p^ = 1 in our model). The effects of 
setting < 1, however, are drastic, and appear for sufficiently strong coupling a, when 
we have observed that oscillations vanish. This is in contrast with the transition to 
the absorbing state found in Ref. |20j for strong coupling. While in their model the 
transition is due to very large amplitudes driving the system into rest (-Pt(l) = 0), in 
our model the system is thrown into an active albeit disordered state, which still has 
Pt(l) 7^ 0, but no oscillations. 

Furthermore, only for non-deterministic excitable elements do we observe 
bistability, with an oscillating and an active (but non-oscillating) phase coexisting. This 
leads to hysteresis cycles, whose sizes can depend on the system size (notably for the 
complete graph) and eventually disappear in random graphs with small enough mean 
connectivity K. 

Although we have restricted ourselves to r = 3, preliminary results suggest that the 
overall scenario is preserved for larger values of r, specially regarding the first transition 
at cTc- The observation of bistability and hysteresis is more difficult for larger values of 
r, owing to stronger finite-size effects. These are similar to those reported by Girvan et 
al.: for finite N and sufficiently strong coupling, the GIG grows in amplitude and nears 
the absorbing state, to which the system is thrown by fluctuations [20]. Distinguishing 
between that type of transition and the bistability reported here is not obvious and 
remains to be studied. 

It is interesting to note that a simple model allows the straightforward application 
of standard techniques of nonlinear dynamics to the study of these phase transitions, 
which could otherwise be difficult to tackle. Note that in the limit = 1, each excitable 
unit in our model has, at the end of its refractory period, a perfect memory of its past r 
time steps. Incorporating this memory in a continuous-time model would require non- 
Markovian dynamics, as recently proposed by Gongalves et al. [30]. Interestingly, their 
model also shows an oscillating phase with reentrance, whose size decreases as memory 
time decreases. In our model, this corresponds to lowering p^, with a similar effect on 
the collective behavior. It remains to be investigated whether bistability also occurs in 
non-Markovian continuous-time models. 

Finally, a note of caution is in order regarding the scaling results of figure [71 The 
fact that a data collapse is obtained employing an upper critical dimension (i^ = 4 by no 
means implies that a lower critical dimension exists. In fact, we are not aware of models 
which exhibit collective oscillations of excitable elements in hypercubic lattices (though 
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they do appear if stimulated with a Poisson drive [IH])- It remains to be investigated 
whether the results of this model hold for small-world networks and other complex 
topologies [2T] . which are extremely appealing for applications in Neuroscience. 
This is currently under investigations (results will be published elsewhere). 

In summary, we have shown that collective oscillations of excitable elements have 
robustness to a certain degree of stochasticity in their intrinsic dynamics. For fixed 
coupling, there is a critical value of p^, below which no oscillations are stable. Fixing 
< 1 and increasing the coupling a, on the other hand, leads to interesting new 
phenomena such as bistability and discontinuous transitions. Taken together, our 
results suggest that even weakly noisy dynamics can qualitatively change the collective 
oscillating behavior. Studies attempting to verify whether these phenomena are observed 
in more detailed excitable networks (e.g. modelled by stochastic differential equations) 
would certainly be welcome. 
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Appendix A. Lyapunov coefficient 



Let u and v be respectively the right and left (adjoint) eigenvector of the Jacobian 
matrix: 

Au =e^^0M, (A.l) 

A^v = e-''^''v , (A.2) 

with both normalized: ({7, u) = 1 = {u, u) (brackets denote the standard complex inner 
product). Let also B{x,y) and C{x,y,z) be multilinear functions proportional to the 
first nonlinear terms of the Taylor expansion of F at a = a^^ , i.e.. 



k,l=l 



XkVi 

NS\ 



(A.3) 
(A.4) 



If we now define 

f = {I - A)-^B{u,'u) , 

s = {e^'U - AY^B{u,u) , 

where / is the t x t identity matrix and u is the conjugate of u, the coefficient li is 
finally given by [12] 



(A.5) 
(A.6) 



V, C{u, u,u)) + 2{v, B{u, r)) + {v,B{u,s)) . (A.7) 
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